Is dislocation flow turbulent in deformed crystals? 
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Intriguing analogies were found between a model of plastic deformation in crystals and turbu- 
lence in fluids. A study of this model provides remarkable explanations of known experiments and 
predicts fractal dislocation pattern formation. Further, the challenges encountered resemble those 
in turbulence, which is exemplified in a comparison with the Ray leigh- Taylor instability. 



From horseshoes and knives to bridges and air crafts, 
mankind has spent five millennia studying how the struc- 
tural properties of metals depend not only on their con- 
stituents, but also how the atoms are arranged and re- 
arranged as metals are cast, hammered, rolled, and bent 
into place. A key part of the physics of this plastic dis- 
tortion is played by the motion of intrinsic line defects 
called dislocations, and how they move and rearrange to 
allow the crystal to change shape. 

Here, we describe the intriguing analogies we found 
between our model of plastic deformation in crystals 
and turbulence in fluids. Studying this model led 
us to remarkable explanations of existing experiments 
and let us predict fractal dislocation pattern formation. 
The challenges we encountered resemble those in turbu- 
lence, which we describe here with a comparison to the 
Ray leigh- Taylor instability. 

For brevity, we offer a minimal problem description 




(a) Continuum Dislocation Dynamics (b) Turbulence 

FIG. 1: Comparison of our continuum 
dislocation dynamics (CDD) with turbulence. 

(a) Dislocation density profile as it evolves from a 
smooth random initial condition. The structures form 
fractal cell wall patterns. Dark regions represent high 
dislocation density. |(b)| Rayleigh- Taylor instability at a 
late time. The fluid (air) with two layers of different 
densities mix under the effect of gravity. The emerging 
flows exhibit complex swirling turbulent patterns. The 
color represents density (red for high, blue for low). 



that ignores many important features of plastic deforma- 
tion of crystals, including yield stress, work hardening, 
dislocation entanglement, and dependence on material 
properties [1]. We focus on the complex cellular struc- 
tures that develop in deformed crystals, which appear to 
be fractal in some experiments [2j. These fractal struc- 
tures are reproduced by our continuum dislocation dy- 
namics (CDD) [3] theory (see Figure [Ta]) . 

Not only do the resulting patterns match the experi- 
mental ones, but the theory also has rich dynamics, akin 
to turbulence. This raises a question: Is the dislocation 
flow turbulent? Here, we focus on exploring this question 
by building analogies to an explicit turbulence example: 
the Rayleigh- Taylor instability. As we describe, our the- 
ory displays similar conceptual and computational chal- 
lenges as does this example, which reassures us that we're 
on firm ground. 

This CDD model [3-5] provides a deterministic expla- 
nation for the emergence of fractal wall patterns jjj 0] 
in mesoscale plasticity. The crystal's state is described 
by the deformation-mediating dislocation density Qij - 
where i denotes the direction of the dislocation lines and 
j their Burgers vectors [l| - and our dynamical evolu- 
tion moves this density with a local velocity Ve, yielding 
a partial differential equation (PDE): 



dtQij 



id m (s n e k VeQkj) = vd i g^ 



(1) 



Here Ve is proportional to the net force on it (overdamped 
motion), coming from the other dislocations and the ex- 
ternal stress. That is, 



Ve 



D 



l 0~mn££mkQkn 



where a is the local stress tensor, the sum of an external 
stress a^ xt and the long-range interactions between dis- 
locations. a\f = J K ijmn (r - r')Q mn (r')dr', with K ijmn 
the function representing the stress at r generated by g 
at r' [1]. The term proportional to v is the regularizing 
quartic diffusion term for the dislocation density (an ar- 
tificial viscosity), which we'll focus on here. (In fact, the 
equation we simulate here is further complicated to con- 
strain the motion of the dislocations to the glide plane 
while minimizing the elastic energy jjjllil.) The details of 
our equations aren't crucial: dislocations move around 
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with velocity V, pushed by external loads and internal 
stresses to lower their energies. Our equation is nonlin- 
ear, and it's exactly this non-linearity that makes our 
theory different from more traditional theories of contin- 
uum plasticity. 

Turbulence is an emergent chaotic flow, typically de- 
scribed by the evolution of the Navier-Stokes equations 
at high Reynolds numbers: 

p(d t v + v-Vv) p\7 2 v + f (2) 
d t p + V-(pv) =0 { } 

where p is the local density of the fluid with velocity v 
under the application of local external force density /. 
The term proportional to p is the fluid viscosity, and p 
is inversely proportional to the Reynolds number. 

Despite this Navier-Stokes equation's enormous suc- 
cess in describing various experiments, there are many 
mathematical and numerical open questions associated 
with its behavior as p — >• 0. In this regime, complex 
scale- invariant patterns of eddies and swirls develop in a 
way that isn't fully understood: turbulence remains one 
of the classic unsolved problems of science. 

How is v related to pi Eq. ([2]) can be written differently 
by dividing the whole equation by p, in this equation 
p/p (in the incompressible case) is called the kinematic 
viscosity (usually denoted as v). Our artificial viscosity 
v in Eq. (pQ) is analogous to this kinematic viscosity. For 
turbulence, p in Eq. ([2]) is given by nature. In contrast, 
our v in Eq. (pQ) is added for numerical stability; it smears 
singular walls to give regularized solutions. This is phys- 
ically justified because the atomic lattice always provides 
a cutoff scale. How do we know that this artificial term 
gives the 'correct' answer (given that there can be many 
different solutions to the same PDE)? Numerical meth- 
ods for shock-admitting PDEs are validated by showing 
that the vanishing grid spacing limit h —> gives the 
same solution as the v — » limit (that is, the viscosity 
solution). We will argue that both our model of plasticity 
and the Navier-Stokes equations do not have convergent 
solutions as p or v — >> 0. We're reassured from the turbu- 
lence analogy and are satisfied with extracting physically 
sensible results from our plasticity theory - viewing it not 
as the theory of plasticity, but as an acceptable theory. 

To solve Eq. (pQ), we implemented a second order cen- 
tral upwind method [6] especially developed and tested 
for conservation laws, such as Eqs. (pQ) and (|2j). The 
method uses a generalized approximate Riemann solver 
which doesn't demand the full knowledge of characteris- 
tics 0. For the simulations of Navier-Stokes dynamics 
(Eq. (|2J)) we use PLUTO [7], a software package built to 
run hydrodynamics and magnetohydrodynamics simula- 
tions, using the Roe approximate Riemann solver. 

Accurately capturing singular flows is a challenge in 
computational fluid dynamics. A classic example of such 
singularities is the sonic boom that happens when an ob- 
ject passes through a compressible fluid (described by a 
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(a) Non-convergence in dislocation dynamics (see Eq. ([TJ) 




(b) N on- convergence in turbulence 
dynamics ( Ray leigh- Taylor instability) 



FIG. 2: Non-convergence exhibited in both 
plasticity and turbulence. As time progresses, the 
curves, which are initially monotonically decreasing, flip 
order and become non- convergent (where the lines cross 
each other). Red arrows show where the convergence is 
lost. 



version of Navier-Stokes) faster than its speed of sound. 
The sonic boom is a sharp jump in density and pres- 
sure, which causes the continuum equations to become 
ill-defined. Our PDEs, depending on gradients of g, be- 
come ill-defined when g develops an infinite gradient at 
a dislocation density jump. 

The numerical methods we use are designed to appro- 
priately solve the so-called Riemann problem: the evo- 
lution of a simple initial condition with a single step in 
the conserved physical quantities. For hyperbolic conser- 
vation laws, exact solutions of the Riemann problem can 
be obtained by decomposing the step into characteristic 
waves. However, in most non linear problems, finding ex- 
act Riemann solutions involves iterative processes th at 
are either slow or (practically) impossible, and thus ap- 
proximate solutions are employed instead. Both methods 



3 



we use are approximate in different ways, but are quali- 
tatively similar. 

These sophisticated methods are designed to handle 
the kind of density jumps seen in sonic booms. In our 
dislocation dynamics, though, we have a more severe sin- 
gularity that forms - a sharp wall of dislocations that 
becomes a ^-function singularity in the dislocation den- 
sity (as v — >• 0). These (^-shocks are naturally present in 
crystals - they describe, for example, the grain bound- 
aries found in poly crystalline metals, which (in the con- 
tinuum limit) form sharp walls of dislocations separating 
dislocation-free crystallites. Unfortunately, the mathe- 
matical and computational understanding of PDEs form- 
ing (5-shocks is relatively primitive; there are only a few 
analytic and numerical studies in one dimension (for an 
example, see [8]). Currently, to our knowledge, there's 
no numerical method especially designed for (5-shock so- 
lutions. Moreover, in a strict mathematical sense, several 
properties of Eq. (pQ) - nonlocality and mixed hyperbolic 
and parabolic features - haven't been proven to permit a 
successful application of shock-resolving numerical meth- 
ods. 

In the simulations presented here, we won't add an 
explicit viscosity (so \i = v = 0); instead, we have an 
effective numerical dissipation [6] that depends on the 
grid spacing h as h n , where n depends on the numerical 
method used. (Eq. (j2j) with fi = is the compressible 
Euler equation. Although we present our simulations as 
small numerical viscosity limits of Navier- Stokes, they 
could be viewed as particular approximate solutions of 
these Euler equations.) 

Figure [T] shows typical emerging structures in simu- 
lations of both the CDD Eq. (pQ) and the Navier-Stokes 
dynamics Eq. (|2j): both are complex, displaying struc- 
tures at many different length scales; sharp, irregular 
walls in the CDD and vortices in Navier-Stokes. Al- 
though it might not be surprising to professionals in 
fluid mechanics that nonlinear PDEs have complex, self- 
similar solutions, it's quite startling to those studying 
plasticity that their theories can contain such complex- 
ity (even though this complexity has been observed in ex- 
periments [!]): traditional plasticity simulations do not 
lead to such structures. 

These rich and exotic solutions demand scrutiny. How 
do we confirm the validity of our solutions? For contin- 
uum PDEs solved on a grid, an important problem that 
needs to be addressed is the effect of the imposed grid. 
Traditionally, it's expected that as the grid becomes finer, 
the solution is likely to be closer to the real continuum so- 
lution. For differential equations that generate singulari- 
ties, one cannot expect simple convergence at the singu- 
lar point! How do we define convergence when singulari- 
ties are expected? For ordinary density-jump shocks like 
sonic booms, mathematicians have defined the concept of 
a weak solution: it's a solution to the integrated version 
of the original equation, bypassing singular derivatives. 



For many problems, researchers have shown that 
adding an artificial viscosity and taking the limit to zero 
yields a weak solution to the problem. For some prob- 
lems, the weak solution is unique, while for others there 
can be several: different numerical methods or types of 
regularizing viscosities can yield different dynamics of the 
singularities. This makes physical sense: if a singular de- 
fect (a dislocation or a grain boundary) is defined on 
an atomic scale, shouldn't the details of how the atoms 
move (ignored in the continuum theory) be important 
for the defect's motion? In the particular case of sonic 
booms, the microscopic physics picks out the viscosity 
solution (given by an appropriate /i — ^ limit), lead- 
ing mathematicians to largely ignore the question of how 
micro-scale physics determines the singularities' motion. 

However, our problems here are more severe than pick- 
ing out a particular weak solution. Neither our disloca- 
tion dynamics nor the Navier-Stokes equation (with very 
high Reynolds number) converge in the continuum limit 
even for gross features, whether we take the grid size to 
zero in the upwind schemes or we take v — ^ (or ji — > 0) 
as a mathematical limit. 

Figure l2al shows a quantitative measure of the our sim- 
ulation's convergence as a function of time, as the grid 
spacing h = 1/N becomes smaller. We measure conver- 
gence using the L2 norm 

where g 2N has been suitably smeared to the resolution 
of g N . (Normally we'd check the difference between the 
current solution and the true answer, but here we don't 
know the true answer.) Here, We study the relaxation of 
a smooth but randomly chosen initial condition - that is, 
a perfect single crystal beaten with mesoscale hammers 
with round heads - as a function of time. We see that 
for short times these distances converge rapidly to zero, 
implying convergence of our solution in the L2 norm. 
However, at around t = 0.02 to 0.2, the solutions begin to 
become increasingly different as the grid spacing h — » 0. 

This worried us at the beginning because it suggests 
that the numerical results might be dependent somehow 
on the artificial finite- difference grid we use to discretize 
the problem, and therefore might not reflect the correct 
continuum physical solution. We checked this by adding 
the aforementioned artificial viscosity v in Eq. (pQ). We 
found that it converges nicely when v is fixed as the grid 
spacing goes to zero. However, this converged solution 
is not unique: it keeps changing as v — >• 0. So, it's our 
fundamental equation of motion (Eq. (pQ)) and not our nu- 
merical method that fails to have a continuum solution. 
This would seem even more worrisome: How do we un- 
derstand a continuum theory whose predictions seem to 
depend on the smallest studied length scale (the atomic 
size)? 



(a) h = 1/128 



(b) h = 1/256 



(c) h = 1/512 



FIG. 3: Continuum dislocation dynamics. Simulation results at t = 1.0 of our CDD Eq. (pQ) at different grid 
sizes (h), starting from a smooth initial condition. We use periodic boundaries in both horizontal and vertical 
directions, and all physical quantities are constant along the perpendicular direction. 




(a) h = 1/128 (b)h = 1/256 (c) h = 1/512 



FIG. 4: The Ray leigh- Taylor instability of the 
Navier- Stokes equation. The Ray leigh- Taylor 
instability is a fluid mixing phenomenon that occurs 
when an interface between two different fluid densities 
is pulled by gravity. These simulation results here are 
the Ray leigh- Taylor instability at t = 4.0, for \i — » 0, at 
different grid spacings (K) with periodic boundaries in 
the horizontal direction and fixed boundaries along the 
vertical. The initial condition has density interface with 
a single mode perturbation in the vertical velocity. The 
system size is (L x ,L y ) = (1.0,2.0). 



It's here that the analogy to turbulence has been cru- 
cial for understanding the physics. It's certainly not ob- 
vious that the limit of strong turbulence \i — >• in Navier- 
Stokes (Eq. ((SJ)) should converge to a limiting flow. Ac- 
tually, our short experience suggests that there is no vis- 
cosity solution for Eq. ([2]). Turbulence has a hierarchy of 
eddies and swirls on all length scales, and as the viscos- 



ity decreases (for fixed initial conditions and loading) not 
only do the small-scale eddies get smaller, but also the 
position of the large-scale eddies at fixed time change as 
the viscosity or grid size is reduced. 

Figure [20 depicts the convergence behavior of a simu- 
lation of the Rayleigh- Taylor instability (as in Fig. l2aj) . 
The instability triggers turbulent flow, and like our dis- 
location simulations, convergence is lost after t ~ 2.0 
as the grid spacing h gets smaller. Our choice of the 
Rayleigh- Taylor instability for comparison is motivated 
by the presence of robust self-similar features (such as the 
"bubbles" and "spikes" in Fig. H] [9]), and by the spatio- 
temporally non-converging features of the initially well- 
defined interface. 

The interface between two fluid densities is analogous 
to our dislocation cell walls. Even though the Rayleigh- 
Taylor instability is different from homogeneous turbu- 
lence in important ways, we also verified that the lat- 
ter shows similar spatio-temporal non-convergence but 
statistical convergence (simulating the Kelvin-Helmholtz 
instability for compressible flow in 2D). The Rayleigh- 
Taylor instability provides the best visualization of the 
analogy between the two phenomena, but this non- 
convergence appears to be more general. Fig. |4] shows 
density profiles at intermediate times for the Rayleigh- 
Taylor instability. The figure shows the formation of 
vortices (swirling patterns), and the simulations look sig- 
nificantly different as the grid spacing decreases. Again, 
this is analogous to the corresponding simulations of our 
dislocation dynamics shown in Fig. [3j where larger cells 
continue to distort and shift as the grid spacing decreases. 

If the simulations aren't convergent, how can we decide 
if the theory is physically relevant and can be trusted 
to interpret experiments? In turbulence, it has long 
been known that, as vortices develop, self-similar pat- 
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FIG. 5: Statistical properties and convergence 

Although the continuum dislocation dynamics (CDD) 
simulations with different grid sizes are 
no n- convergent (Fig. [2aj), the statistical properties are 
the same. The dislocation density correlation function 
is plotted here for different simulation sizes at the same 
time, exhibiting consistent power laws. 



terns arise in the flow and exhibit power laws in the 
energy spectrum and in the velocities' correlation func- 
tions [10]. A successful simulation of fully developed tur- 
bulence isn't judged by whether the flow duplicates an 
exact solution of Navier-Stokes! Turbulence simulations 
study these power laws, comparing them to analytical 
predictions and experimental measurements. 

Our primary theoretical focus in our plasticity study [3] 
has been to analyze power-law correlation functions for 
the dislocation density, plastic distortion tensor, and lo- 
cal crystalline orientation. As Fig. [5] shows, like tur- 
bulence simulations, these statistical properties seem to 
converge nicely in the continuum limit. 

It's worth noting that, in both cases, non-convergence 
emerges when small scale features appear on the wall (see 
Fig. [6]) or the interface (see Fig. [7j): 

In the case of our simulations of plastic flow (Eq.[T]), 
starting from smooth initial density profiles, finite time 
singularities develop in the form of J-shocks. The ex- 
istence of finite-time singularities was shown in a ID 
variant of these equations, which is associated with the 
Burgers equation [11]. Figure [6] shows how this effect 
occurs by considering the two-norm differences (the in- 
tegrand in space of the norm discussed earlier). At 
t = 0.04 (see Fig. [2aj) when the N = 512 curve starts to 
cross all the other curves, singular features start to ap- 
pear around a wall (Fig. l6b|) . Although the boundaries 
are no n- convergent when specific locations and times are 
considered (Fig. [2aJ), the statistical properties and asso- 
ciated self-similarity (Fig. [5]) are convergent. 

In the case of Navier-Stokes simulations (see Fig. [7j), 



the existence of finite-time singularities is a topic of active 
research: even though local-in-time analytic solutions are 
easily shown to exist, global-in-time analytic solutions 
can be proven to exist only for special cases, such as in 
the 2D incompressible genuine Euler equation [Hj]. I n 
3D, the mechanism of vortex stretching is conjectured to 
lead to finite-time singularities [13], even though there 
are still crucial open questions. Despite its complexity, 
turbulence can be concretely studied in special cases. For 
our example of the Ray leigh- Taylor instability, the two- 
fluid interface gets distorted and "bubbles" form (Fig.H]); 
over time, the bubbles exhibit emergent self-similar char- 
acteristics [9], showing statistical convergence. However, 
there is no spatio-temporal convergence, because the in- 
terface develops complex, turbulent features as the grid 
becomes finer (see Figs. I2blandl7j). 

Sometimes science seems to be fragmented, with in- 
dependent fields whose vocabularies, toolkits, and even 
philosophies almost completely separate. But many valu- 
able insights and advances arise when ideas from one field 
are linked to another. Computational science is provid- 
ing a new source of these links, by tying together fields 
that can fruitfully share numerical methods. 

Our use of well-established numerical methods from 
the fluids community made it both natural and easy to 
utilize their analytical methods for judging the validity 
of our simulations and interpreting their results. 
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(a) h = 1/256 



(b) h = 1/512 




(c) h = 1/256 



(d) h = 1/512 



t = 0.04 



t = 1.00 



FIG. 6: Non-convergence and singularity for CDD For [(a)] and [(b)) £ = 0.04; for [(c)] and p)| t = 1.00. The 

two-norm difference between h and h/2 are plotted. At short time, t = 0.04, the differences are small: |(a)| is empty 
and |(b) | nearly so. At later times, t = 1.00, the two-norm difference becomes significant esepcially where the walls 
are forming (see Fig. [3] for wall locations). 




(a) h = 1/64 



(c) h = 1/256 




(d) h = 1/64 



(e) /i = 1/128 
t = 3.0 




(f) h = 1/256 



FIG. 7: Non-convergence and vortices for Navier-Stokes For (a) , |(b) , and |(c)| t = 2.0; for (d) , [(e), and |(f)| 
t = 3.0, corresponding to the two red arrows in Figure l2bl The two -norm difference between h and h/2 is plotted. 
At t = 2.0, small structures start to become strong, as in |(c)[ and the same is true at t = 3.0 for h = 1/128, as in [(e) 
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